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Abstract 

The microscopic basis of Newton's law of cooling and its modification when the differ- 
ence in temperature between the system and the surroundings is very large is discussed. 
When the system of interest is interacting with a small bath, the effect of the dynamical 
evolution of the bath variables is important to find out its dynamical feedback on the 
system. As in the usual system-bath approach, however, the bath is finally considered to 
be in thermal equilibrium and thereby provides an effective generalization of the Born- 
Markov master equation. It is shown that the cooling at early time is faster than that 
predicted by Newton's law due to the dynamical feedback of the bath. 
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I. Introduction 



Classically the phenomenon of cooling of a bulk body may be considered as a process where 
the flow of heat from the hotter body to a colder environment is governed by the Newton's law 
of cooling, namely, 

- d -W = l(T(t)-T R ), (1) 

with T(t) be the instantaneous temperature of the body, Tr be the temperature of the en- 
vironment and 7 be the characteristic decay constant, respectively. The solution of Eq.(l) 
reads 

r(t)=T e-T* + r fl (l-e-^), (2) 

where T is the initial temperature of the body at t — such that T R < T . The classical 
rate equation does not contain any non-linear term and it is valid if the difference between the 
temperature of the system and environment is small. 

On the other hand, the cooling of the neutral atoms is generally performed by various 
coherent laser-cooling techniques, which primarily concern with the reduction of the kinetic 
motion of the center-of-mass of the trapped atoms [1,2]. However, if we consider the coherent 
cooling of the molecules [3], then the contribution to the thermal energy arises not only from the 
motion of their centre of mass, but also from their rotational or vibrational motion. Thus for an 
ensemble of trapped molecules the cooling requires the ceasing of the momentum associated with 
all possible degrees of freedom, although a priori it is difficult to ascertain which of the degree 
of freedom contributes most significantly. If we neglect the translational and rotational motion 
of the trapped molecules, the dissipation of temperature associated with the vibrational degrees 
of freedom may be formulated quantitatively within the framework of the standard dissipation 
theory of the damped harmonic oscillator. To formulate the theoretical basis of the vibrational 
cooling, we assume that a molecule with one or few modes of vibration as our system of 
interest which dissipates its energy into the large number of other modes acting effectively acts 
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as the reservoir. In the density matrix formalism of the system-reservoir composite model, the 
nontrivial weak coupling between the system with the reservoir effectively induces a damping in 
the system. Consequently, we obtain the Born-Markov master equation of the reduced density 
matrix p, where the reservoir oscillators are completely eliminated in terms of the system 
variables [4,5] as 



^ = -iu [ata, p] - 7(1 + n R (T R ))(a) ap - 2apa t + pa)a) - jn R (T R )(aatp - 2a) pa + paa 1 "), (3) 

where the frequency of the system oscillator is u with [a, a"!"] = 1 and the average thermal 

excitation number of the bath is n R (T R ) = — , where T R is the temperature of the 

- 1 

reservoir. Thus, it follows from the master equation that the time evolution of average photon 
number of the system mode, dn ^ = Tr[a)a^\ with frequency loq is described by the rate 
equation 

dn s {t ) 
dt 

and its solution is given by 



= -y(n s (t)-n R (T R )), (4) 



n s (t) = n R (T R ) + (n s (0) - n R (T J? ))e^ t . (5) 

Thus, the vibrational cooling corresponds to the feeding of the thermal photon from the system 
to the reservoir until the system photon number equilibrates with that of the reservoir. This 
process is known as the thermalization and it is evident from Eq.(5) that, similar to the classical 
cooling, it occurs after an infinitely large time, namely, 77.5(00) = n R {T R ). We can associate the 
instantaneous average photon number of the system mode with an effective temperature T(t) 
by using the relation 

ns(T(t)) = , (6) 

which in the high temperature limit gives, n s (T(t)) ~ In the same limit, the initial system 

photon number and the reservoir photon number are given by ns(To) |^ and n R (T R ) as 



respectively; putting these values in Eq.(5) we recover Eq.(2). Therefore, the classical Newton's 
cooling law in Eqs.(l) and (2) can be obtained from the high temperature limiting situation 
of the Born-Markov master equation in Eqs.(3) and (5). This is valid when the difference in 
average energy per mode between the system and the reservoir is small, whereby the reservoir 
does not change with time due to the acceptance of energy from the system of interest. 

In the crossroad of various approaches and applications of the system-reservoir composite 
formalism, the necessity of the finite bandwidth of the reservoir [6-12] was envisaged right 
from the beginning which leads to the possible modification of Eq.(3). The recent experiments 
involving ultra-fast time scale [13-15], correlated emission laser (CEL) pulse with adjustable 
memory time [16], experiments on cavity electrodynamics [17] etc, have significant impact in 
the understanding of the models beyond the Born-Markov approximation. However, in all 
previous studies, the assumption that works at more subtle level is that the photon absorbed 
by the reservoir from the system cannot bring any dynamical change with it because of the 
small difference of the energy between the system and average energy of each degree of freedom 
of the reservoir. It is therefore of natural interest how the situation changes if we consider that 
the energy of the system is large enough in comparison to that of each mode of the reservoir. 
In this paper we shall show that, due to the large difference of average energy between the 
system and each reservoir mode, the flow of the thermal photon of substantial energy from 
the system to the reservoir effectively leads to the dynamical evolution of the reservoir. Our 
primary objective is to discuss the vibrational cooling of a system with large energy content by 
incorporating the aforesaid dynamical evolution of the reservoir. 

The remaining sections of the paper are organized as follows. In Sec. II we have developed 
a formalism to incorporate the evolution of the reservoir and show how it effectively generalizes 
the Born-Markov master equation beyond leading order of the decay constant. We apply our 
formalism to the damped harmonic oscillator in Sec. Ill and show that how it affects the 
thermalization time to make it short. We conclude by highlighting the essence of the paper 
and discussing its outlook. 
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II. Formalism 

The Hamiltonian of a system interacting with the reservoir is given by 

H = H S + H R + V = H + V, (7) 

where H$, Hr represent the Hamiltonians of the system and reservoir, respectively and V is 
the interaction between them. Let n{t) be the joint density matrix of the composite system in 
the Interaction Picture (IP); the corresponding evolution equation is 

^ = ><*),*>]■ (8) 
The solution of the equation is given by 

K{t) = «(0) - ifdt>[V(t>),K(0)} + (f ) 2 }dt'fdt"[V(t'), [V (?'),*(?')}}■ (9) 



We consider the interaction Hamiltonian in the IP is of the following form: 

v{t)=hZQi(t)Fi(t), (10) 

i 

where Qi{t) and Fi(t) be the system and reservoir operators, respectively, in IP. Using the 
factorization ansatz, namely, k(0) « s(0)/(0) and n(t") s(t")f(t"), respectively, and by 
noting the fact that Tr R n(t) = s{t), the trace over the reservoir mode in Eq.(9) yields 

s(t) = s(0) - ifdt f E(Fi(t'))RlQi(t'),sm ~ 

i 

j dt' j dt'' J2{Tr R [F i (t')F j (t' / )f(t'')][Qi(t / )Q j (t / ')s(t'') - Q j (t")s(t")Q i (t')]- 

o o i,j 

TrniF^F^fit^m^sinQAn ~ ^'OW)^')]}, (11) 

where {...)r = TV\r[.../(0)] is the average of the reservoir operators. Taking the derivative of 
Eq.(ll) with respect to t' we obtain (we redefine t" by t' and t' by t), 

^ = -£/ dt'iUQ^t') - QX*)Qi]TrRm)FiV)fV)]- 

[Q i s{t')Q j - sWQ^TrnlF^F^fimexp^t + cuff), (12) 
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where uf be the characteristic frequencies of the system and the term linear in reservoir operator 
vanishes by the symmetry argument. The system oscillator in the IP in V(t) can be expressed 
in the Schrodinger picture (SP) by using the standard prescription, 

Q t (t) = e m)H S tQ. e -{im{H S yt 

= Q t er^\ (13) 

Now, replacing t' by t — r in Eq.(12) and assuming the Born-Markov approximation [4], namely, 
s(t — t) ~ s(t) for large t, the generalized master equation of the reduced density operator in 
Schrodinger picture, S, is obtained as 

f = -i[H s ,S] -EmQjS-QjSQ^W^t]- 

[Q.SQ, - SQ^W-MW"? + uf), (14) 

where 

W±[t}=jdre^Tr R [F l (t)F J (t-T)f(t-r)], (15a) 

o 

W'M =Jdre^Tr R [F 3 (t - r)F t (t)f(t - r)], (15b) 

which are to be calculated in different situations. 

To include the evolution of the reservoir in this scenario, using Trs^{t) = f(t) along with 
Eq.(lO), the trace of Eq.(9) over the system yields 

f(t -r) = /(0) - i E 7 <fti<Q*(*i))s[**(*i), /(0)]- 

k 

E Tdh /dt 2 [(F,(ti)F m (t 2 )/(0) - F m (t 2 )/(0)^(ti))(QK*i)<5m(t2)>5- 

l,m 

(F,(tO/(0)F m (t 2 ) - /(0)F m (i 2 )F I (*i))(g m (* 2 )Qi(*i)>s] + -> (16) 

where (...)s = Trs[...s(0)] represents the average of the system operator which depends on 
the initial population distribution of the system. In deriving Eq.(16), unlike previous case, we 
use the ansatz n(t") k(0) to terminate the series beyond the second order in interaction [4]. 
Plucking Eq.(16) back into Eq.(15a), we find 
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W$[t]=W*°[ t ] + W$[t] + ..., (17) 

where = / dre 1 ^ 7 '{Fi(t)Fj(t — T))R is the usual lowest order reservoir correlator [4,5]. In 

o 

Eq.(17) the term next to lowest order arises due to the correlation among the system oscillators, 
and it is given by 

W$[t] = -iEldre^f^TdhiQkih^siiF^Fjit-^Fkit^R- (F^F^F^ - r)) R }- 

k 

E J dre iul i T J dt l fdt 2 

l,mO 

[((F^F^t - r)F,(*i)F ro (t 2 )) fl - (FfrjFiWFjit - r)F m (t 2 )) R )(Q l (t 1 )Q m (t 2 )) s + 
( ( F m (t 2 )Fi (ti)Fi(t) Fj (t - t)} r - (F^F^F^t - r)F,(ti)) fl )(g ro (t 2 )Qi(*i))s]- (18) 

It is customary to write the reservoir Hamiltonian in the following form: 

H R = Ef^Mb k + l), (19) 
k 

where Uk is the frequency of the reservoir modes. The time-dependent reservoir operator in the 
IP appearing in Eq.(18) can be expressed in the SP as 



Fi(t) = E K p e {i/h)HRt b p e ( ' i/h)HRt 
v 

= E Kp&pe-*"**, (20a) 
v 

F 2 (t) =Y,K q e^ h)HRt bU- i /V H * t 
i 

= E KMe^, (20b) 

where lo\ (I = p, q) is the angular frequency of the reservoir oscillators mode and Hi the coupling 
constants, respectively. From Eq.(18) we now proceed to evaluate the spectral density function 
W&[t] for % = 1 and j = 2, 

W+[t] = ~ E I dre^T dh Jdt 2 

LmO 



[((F 1 (t)F 2 (t - r)F,(*i)F m (t 2 )) fl - (F l (t 1 )F 1 (t)F 2 (t - r)F m {t 2 )) R ){Q l (t 1 )Q m (t 2 )) s + 

((F m (t2)F l (t 1 )F 1 (t)F 2 (t - t)) r - (F m (t2)F 1 (t)F 2 (t - ^F^r^QMQ^s], (21) 

where, once again, the term linear in the system operators is dropped by the symmetry argu- 
ment. Throughout the treatment we assume that the reservoir is in a thermal distribution and 
thus only the diagonal terms will survive. Substituting Eqs.(13), (20a) and (20b) in Eq.(21), 
we obtain (for details, see Appendix) 

W&[t] = E |K r | 2 |K*| 2 e i(Wr ~ Ws) * f dTe i( - u °- Ur > r Y dhe^-^* 1 

r,s 

/ dt 2 e^~^{{2 + n R {u r ) + n R (u s ,T R ))(Q 1 Q 2 ) s - (n fl (w r ,T fl ) + n fl (w a ,T fl ))(g 2 gi) s ), (22) 
o 

where n R (u s , T R ) is the average thermal photon number of the reservoir and we have considered 
the system frequency to be ujf = —uj 2 = ujq for convenience without losing generality. Thus we 
note that in Eq.(22), the evolution of the reservoir arises due to the correlation functions of the 
system operators. Finally converting the sum over modes into the frequency space integrals, 
we find 

Wt 2 \t\ = J du r D{uj r )\K{uj r )\ 2 J du s D(Lu s )\K( y uj s )\ 2 e i ( UJr ~ UJs ^ t J dre'^ ^ QJr ' T 



T^ 1 e^'-^4 1 ^2e^-^ ) * 2 ((2 + n^(cu r ,T^) + n fi (a; s ,T i? ))(Q 1 Q 2 ) 5 - 
o o 

(n R (u r ,T R ) +fi R (co s ,T R )){Q 2 Q 1 )s), (23) 

where D(u r ) and D(u s ) be the density of states respectively. Proceeding in the similar way we 
can show that W^[t] = W 21 [t\. The time development of the reservoir for any simple system 
may be calculated from Eq.(23). 
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III. Application to simple harmonic oscillator 



The vibrational cooling may be modeled by a harmonic oscillator interacting with large 
number of the reservoir modes resulting damping. The free Hamiltonian and interaction term 
of such composite system are given by 

H s = Tiu (cJ a + \) (24a) 

V = hE(K k rfb k + K* k ab{), (24b) 

k 

respectively, with the generic reservoir Hamiltonian given by Eq.(19). In Eqs.(24a) and (24b), 
the system operators are in the Schrodinger picture i.e, Q\ = a) and Q2 = cl- Taking f%s(t) = 
(cra(t))s to be the average photon number of the system in time t and the upper limits of 
integration as t, t — r,t 1 — > 00, Eq.(23) becomes 

W+[t] = 2n 2 D 2 (cu )\K(uj )\\n s (t) - n R (T R ))t, (25) 

the difference between the the instantaneous average excitation number of the system and the 
thermal average photon number of the reservoir. In deriving Eq.(25) we have neglected the 
principal parts which correspond to a small Lamb shift. Substituting Eq.(25) in Eq.(17) (with 
i = 1 and j = 2) and plucking back the resulting equation in Eqs.(15a) and (15b), we obtain a 
generalized Born-Markov master equation of the damped harmonic oscillator, 

f t = -i[H s , S] - ^[rfaS - 2aSa) + Sa^a] - ^[aa^S - 2a^Sa + Saa*}, (26) 

where 

7 i(t) = 7(1 + n R {T R )) + 7 2 (r*s(t) - n R (T R ))t, (27a) 

72W = in R {T R ) + l\n s (t) - n R (T R ))t, (27b) 

with 7 = 27r|K(a;o)| 2 -D(^'o) the decay constant. Thus, a linear time-dependent term appearing 
beyond the leading order of the decay constant becomes important if \{ns{t) — n R (T R ))\ » 0. 
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If the initial average energy of the system is much more than the thermal average excitation 
of the bath, then the time-dependent decay rate 7i(t) and 72 (t) appreciably affect the decay 
dynamics in an early time. 

To address the notion of thermalization in our scenario, we need to calculate the time 
evolution of n s (t) (= (a^a^t^s) from Eq.(26), which is governed by the rate equation 

-^f^ = 7(M*) - n R (T R ))(l + jt) (28) 
and its solution reads 

ns(t) = n R {T R ) + (n s (0) - n R (T R )) e -^ 1+ ^ . (29) 

Comparing Eq.(5) with Eq.(29), we note that in the later case the decay rate is time- dependent 
and thermalization becomes faster due to the incorporation of the dynamical evolution of the 
reservoir. 

Here, we note that in the high temperature approximation, Eq. (28) leads to the modified 
Newton's law of cooling as 

_*m = 7(T(t) _ Tr) + 7 2 (T(t) _ Tr% (30) 

where the term beyond the leading order of the decay constant becomes significant if T >> T R . 
The Fig. shows the comparison of the plots of Eq.(2) with Eq.(30), where we note that the 
thermalization occurs at a faster rate. Cooling in early time is much faster than predicted by 
Newton's law. As a first order correction the theory is valid upto t < 7 _1 ; for a time longer 
than t > 7 _1 , Newton's exponential law should be considered for thermalization. Cooling-time 
for reaching from 2000°C to 800°C is almost two-third in the modified dynamics of that in the 
Newton's cooling law. In comparison to Newton's law, where the time required to bring the 
temperature of the system above the reservoir-temperature to its half, i.e., — T^r = ^ 
half-thermalization time is, ti/ 2 = ^r, whereas for the modified case the half-thermalization 
time is given by ti/ 2 = ^ /1+2 ^ n2 ~ 1 -. Therefore, in the modified case also ti/ 2 is independent of 
(T(0)-T R ). 
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[Fig.]: The plot of Newton's cooling law is given by Eq.(2) (Newton, solid line) and its 
modified form given by Eq.(30) (Modified, dashed line) with time(in 7 _1 ) for initial and 
final temperatures(units arbitrary), T = 2000° C and T R = 200°C. 

From Eq.(26), the corresponding master equation of the diagonal elements of the density 
matrix can be given by the loss-gain equation of population for the step-ladder model of a 
harmonic oscillator as 



dSjjit) 
dt 



= £ [w{j\i)s n {t)-wmsm- 



(31) 



j=i±i 



Here the transition rates connecting only the neighbouring levels are specifically given by 
W(i + l\i) = (i + 1) 7 [1 + {n R + (n s (t) - n R )jt}] 

W(i\i + 1) = n[n R + (n s (t) - n R )jt], (32) 

which means that the (i + l)-th to z-th state transition rate is time dependent, and vice- versa. 
On the top of the usual temperature-dependent rate, it depends on yt and on the difference 
in temperature of the reservoir from that of the instantaneous temperature of the system, 
T s (t) — T R (or equivalently (n s (t) — n R )). Usually we consider < 1 for the first-order 
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perturbative effect. A direct consequence of the loss-gain-type master equation shows that the 
bath- induced forward and backward rates are modified by a factor of n R + {n s {t)—n R )^t instead 
of n R . As the modification arises due to the system-induced dynamics of the reservoir, which is 
considered as a first-order effect, we can safely assume that the rate is primarily governed by the 
factor n R and therefore we can consider {n s {t) — n R )^t < n R . This amounts to the fact that the 
initial system temperature should not be arbitrarily high compared to the temperature of the 
reservoir. Otherwise, a strong non-equilibrium evolution of the bath will produce a nonlinear 
coupled dynamics of system and bath variables, which is immensely difficult problem to tackle 
to provide any tangible physical result. 

IV. Conclusion 

In this paper we have developed the quantum theory of cooling of a system with large 
energy content when the reservoir has also a dynamical evolution rather than thermal equilib- 
rium. It is explicitly demonstrated for such system that the thermal equilibrium is attained 
much faster than in comparison to the case of exponential decay when the reservoir is at equi- 
librium. Our study reveals that, the larger is the initial photon number content of the system, 
the faster is the rate of cooling. Possible modification of the Newton's classical law of cooling 
beyond the leading order of the decay constant is pointed out and it is shown that the higher 
order term becomes significant if the difference between the average energy per mode of the 
system and the reservoir is considerable. The analysis is strictly valid for a very short time, 
t < 7 -1 , and initial temperature or average energy of the system is not too much higher com- 
pared to the reservoir as the modification in the theory arises due to first order perturbation 
effect. We have also considered a repeated neglect of off-diagonal terms corresponding to bath 
degrees of freedom arising from the dynamics where only the diagonal elements of the bath 
density are modified in time. A faster thermalization or faster cooling is qualitatively under- 
standable as the bath is interacting with the system more actively instead of passively waiting 
in its equilibrium distribution in the usual approach. 

To treat a finite size [18,19] of the bath, one immediate choice is to restrict the number 
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of modes in the bath. This is equivalent to a pronounced recurrence of population of the 
vibrational states due to the back and forth exchange of energy between the modes of the system 
and bath. However, in our approach we have assumed the fact that the system experiences a 
feedback due to the dynamical evolution of the bath, but ultimately the bath is assumed to be 
in thermal equilibrium. We have calculated the two-point and four-point correlation functions 
of the bath variables by averaging over the thermal equilibrium distribution instead of a non- 
equilibrium distribution of the bath. The population decay is non-exponential due to this, 
which has simple dependence on the difference between the average energy of each mode from 
the time-dependent state of the system to the bath at equilibrium. A more systematic approach 
to treat finite size of the bath is under consideration which will be published elsewhere. 

In the midst of several currently interesting coherent cooling mechanism of atoms and 
molecules induced by laser, this incoherent mechanism of vibrational cooling may find it worth- 
while because of the huge availability of nano-materials [20] which can support a large number 
of degrees of freedom. It can effectively act as a bath as well as a finite quantum system to 
reciprocates energy with a subsystem of interest which is composed of a single or a few modes of 
vibration. Other associated aspects of the system-reservoir formalism require careful scrutiny 
in the light of the dynamical evolution of the reservoir proposed here. 
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Appendix 

In this appendix we shall derive the four-point reservoir correlation functions appearing in 
Eq.(21). The four-point reservoir correlators can be expressed in terms of two-point correlators 
by using the identity, 

(O a O b O c O d ) R = (O a O b ) R (O c O d ) R + (O a O c ) R (O b O d ) R + (O a O d ) R (O b O c ) R . (A.l) 

Using Eq.(A.l) let us calculate the term in the square bracket of Eq.(21) with % — 1, j — 2 and 
/, m runs from 1 to 2, 

E i , m [(F 1 (t)F 2 (t')FKti)F m (t 2 )) i? -(FKti)F 1 (t)F 2 (t')F m (t 2 )) i? ](Q i g m ) 5 

+[(F m (t 2 )F l (t 1 )F 1 (t)F 2 (t')) R - (F m (t 2 )F 1 (t)F 2 (t')F l (t l )) R ](Q m Q l ) s = 

[(F 1 (t 1 )F 2 (t'))i?((F 2 (t 2 )F 1 (t)) i? - (F l (t)F 2 (t 2 )) R ) + 

(F 1 (t)F 2 (t 2 )) R ((F 2 (t')F 1 (t 1 )) R - (F 1 (t 1 )F 2 (t')) R )](Q 1 Q 2 ) s 

l(F 2 (t' )F 1 (t 1 )) R ((F 1 (t)F 2 (t 2 )) R - (F 2 (t 2 )F 1 (t)) R )+ 

{F 2 {t 2 )F 1 (t)) R ({F 1 {t 1 )F 2 (lf)) R - (F 2 (t')F 1 (t 1 )) iJ )](g 2 g 1 ) 5 . (A.2) 

To obtain the above equation we have neglected the off-diagonal terms, since the reservoir 
assumed to be a thermal one. Substituting Eqs.(13), (20a) and (20b) in Eq.(A.2), the straight- 
forward simplification yields 

E i , m [(F 1 (t)F 2 (t')F(ti)F m (t 2 )) i? -(F(ti)F 1 (t)F 2 (t')F m (t 2 )) i? ](Q i g m ) 5 
+ [(F m (t 2 )F(ti)F 1 (t)F 2 (t / ))i? - (F m (t 2 )F 1 (t)F 2 (t')F(ti))i ? ](Q m ^)5 = 

_ ^ Kir ^^^2 e i(u) r -u) s )t e -i(u) r ti-u) s t2) e -iu) r T 
r,s 

((2 + n R (u r ) + n R (u s )){QiQ 2 ) s - (n R (uj r ) + n R (uj s )){Q 2 Qi) s ), (A. 3) 

where t' is replaced by t — r. Finally substituting Eq.(A.3) in Eq.(21) we obtain Eq.(22). 
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